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Abstract. We investigate the clustering of halos in cosmological models starting with general 
local-type non-Gaussian primordial fluctuations. We employ multiple Gaussian fields and add 
local-type non-Gaussian corrections at arbitrary order to cover a class of models described 
by frequently-discussed f n \, g n \ and r n i parameterization. We derive a general formula for 
the halo power spectrum based on the peak-background split formalism. The resultant 
spectrum is characterized by only two parameters responsible for the scale-dependent bias 
at large scale arising from the primordial non-Gaussianities in addition to the Gaussian bias 
factor. We introduce a new inequality for testing non-Gaussianities originating from multi 
fields, which is directly accessible from the observed power spectrum. We show that this 
inequality is a generalization of the Suyama-Yamaguchi inequality between f n \ and r n i to 
the primordial non-Gaussianities at arbitrary order. We also show that the amplitude of the 
scale-dependent bias is useful to distinguish the simplest quadratic non-Gaussianities (i.e., 
/ n i-type) from higher-order ones (g n \ and higher), if one measures it from multiple species of 
galaxies or clusters of galaxies. We discuss the validity and limitations of our analytic results 
by comparison with numerical simulations in an accompanying paper. 
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1 Introduction 



The large-scale structure (LSS) of the universe as well as the temperature fluctuations of 
the cosmic microwave background (CMB) radiation are well understood in the standard 
framework of cosmological model starting from tiny almost Gaussian fluctuations. However, 
the generation mechanism of the primordial fluctuations, which seed both the CMB and the 
LSS, is yet to be understood. Because we usually consider that the primordial perturbations 
are created during inflation, their statistical properties provide us valuable information about 
the inflationary physics. Although we have no clear evidence of the violation of the Gaussian 
assumption of our initial condition to date, the signature of the primordial non-Gaussianities, 
if they are detected, can be a smoking gun to test various models of the primordial fluctuations 
in the coming era of precision cosmology. 

Some proposed models of the early universe predict that large local-type non-Gaussianities 
are produced when nonlinear dynamics is important on super-horizon scales [1]. In these 
models, the Bardeen's curvature potential <I> in the matter dominant era is written as a local 
function of a Gaussian field eft, and is usually expanded into a Taylor series as 

$(x) = 0(x) + f n i [</> 2 (x) - (0 2 )] + g nl 3 (x) + h nl [^> 4 (x) - (0 4 )] + i nl 5 (x) + • • • . (1.1) 

The coefficient f n \ in the quadratic term is one of a key parameter to distinguish the models 1 
[2]. This parameter determines the amplitude of the primordial bispectrum, the lowest-order 
statistic that captures the deviation from Gaussian statistics. In this model, the bispectrum 
of the curvature perturbation £ = (3/5)$ is written as 

B c (h,k 2 ,k 3 ) = |/ n i[i , c(fei)P f (fc 2 ) + (2perm B .)] I (1-2) 

where P^ denotes the power spectrum of £ and (2 perms.) stands for two more terms given 
as permutations of wavenumbers, k\, k 2 and k 3 . The current constraints on the parameter 
/ n i from the CMB experiments are very tight, and any deviation from Gaussianity of or- 
der ^ 0.1% level in the curvature perturbations has already been ruled out (e.g., [3]). It 
is expected that the constraint will be much tighter after some on-going/upcoming large 
observational projects such as Planck [4]. 

However, the value of / n i does not determine all the statistical properties of the initial 
condition of our universe, even if one focuses only on the local-type non-Gaussian primordial 
fluctuations. The next order statistic, the trispectrum, is a natural extension of the discussion, 
and it indeed gives us more hints about the inflationary physics. The cubic-order correction 
to the curvature potential is a natural source of this statistic [5, 6] (<? n l _ type; see Eq. 1.1). 
This gives the trispectrum of the form 

54 

T c (k 1 ,k 2 ,k 3 ,h) = ^<?ni [P C (h)P c (k 2 )P c (k 3 ) + (3 perms.)] . (1.3) 

The trispectrum can also be generated through the quadratic corrections [7]. This type of 
trispectrum has a different scale dependence from Eq. (1.3) and is written as 

T c (k 1 ,k 2 ,k 3 ,k 4 ) = r n i [P c (fci)P f (fc2)Pf(|ki +k 3 |) + (11 perms.)] , (1.4) 

1 Note that we adopt the normalization of / n i widely used in studies of the CMB, instead of the convention 
sometimes employed in the LSS works, where this parameter is also defined by Eq. (1.1), but the normalization 
of the curvature potential is done at the present by linear extrapolation. 
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where the parameter r n i controls the amplitude of the spectrum. In single-field models 
described by Eq. (1.1), r n i can be expressed by f n \ as r n i = (36/25)/^. In general, however, 
when one employs two or more Gaussian components in the primordial fluctuations, r n i can 
be larger, and the following inequality is shown to be always valid ([8]; Suyama-Yamaguchi 
inequality) : 

-nl > |/ n 2 !, (1.5) 

where f n \ is defined not by Eq. (1.1) but by Eq. (1.2) as the amplitude of the bispectrum of 
the curvature perturbations. See [9-11] for further discussions on the loop corrections and the 
universality of this inequality, and [12] for an observational attempt to test it from existing 
observations. Thanks to this inequality, one can distinguish multi-field inflationary models 
from single-field ones by measuring f n \ and r n i simultaneously. The statistical properties of 
the primordial perturbations are uniquely determined by three parameters, f n \, g n \ and r n i, 
up to the order of trispectrum as long as one is interested in local-type non-Gaussianities 2 , 
and a number of models can be sorted by these parameters [9]. 

The LSS that we observe through the clustering of galaxies or clusters of galaxies is 
a highly non-Gaussian field due to the nonlinear gravitational evolution (e.g., [15] for a 
review), redshift-space distortions caused by the peculiar velocity field of the galaxies [16] 
and the galaxy bias with respect to the underlying matter density field [17]. These acquired 
non-Gaussianities are thought to prevent us to measure only the primordial non-Gaussianity 
from the LSS unlike the CMB, whose temperature fluctuations are still in the linear stage. 
It is thus not straightforward to extract the non-Gaussian signature of purely primordial 
origin from the LSS, and a number of studies have been done to model the gravitationally 
generated non-Gaussianities [18-22]. In these studies, the authors mainly focus on the galaxy 
bispectrum, simply because it is the lowest-order statistic that captures the non-Gaussian 
signature (but see [20] for a discussion on the trispectrum). 

However, one observable feature in the LSS, which was realized only recently, attracts 
a great attention [23]. That is the scale- dependent bias in the halo power spectrum on large 
scales (i.e., ^ Gpc). This is an interesting example where the bias, which usually prevents 
us from extracting cosmological information, opens a new window to see the primordial 
bispectrum but in the power spectrum. This new feature has been extensively studied both 
analytically [24-31] and numerically [32-34], and gives stringent constraints on f n \ [35-37], 
which are already competitive to those from the CMB. Also, the scale-dependent bias affects 
the halo bispectrum similarly, and thus there is a possibility to put even stronger constraints 
on / n i from the LSS by combining the power spectrum with the bispectrum [34, 38-41]. 

A similar scale-dependent halo bias has been reported in the presence of the primor- 
dial trispectrum (see [42-45] for g n i-type and [45-47] for r n i-type). Furthermore, the scale 
dependence in the halo bias is found to be present in other models of the primordial non- 
Gaussianities, such as the scale-dependent f n \ model [48, 49], ungaussiton model [50], and 
non-local non-Gaussian models [51-53]. Thus the interpretation is not straightforward even 
when one detects a clear evidence of the scale-dependent bias from observed clustering of 
galaxies or clusters of galaxies. 

Our aim is to develop a statistical methodology to distinguish these models. Moreover, 
we would like to generalize the argument as model-independent as possible. In this paper, as 
a first attempt, we focus on the scale-independent local- type non-Gaussianities and discuss 

2 Note, however, that these parameters can be scale dependent in general (see, e.g., [13, 14]). 
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the halo clustering. In particular, we employ multiple Gaussian fields that seed the curvature 
perturbations in order to control f n \ and r n i independently. We also keep all the higher-order 
corrections in the curvature perturbations such as g n \ and h n \ to see how generic our result 
is. Our model is a simple generalization of the / n i, g n \ and r n i parameterization. We will 
show that this class of primordial non-Gaussianities result in a universal formula for the 
halo power spectrum. We introduce a new parameterization of the non-Gaussian signature 
in the scale-dependent bias and propose two tests to distinguish different models based on 
the analytic results. The validity and the limitations of our analytical results are tested by 
confronting with a large set of cosmological A-body simulations, and will be presented in a 
separate paper (hereafter paper II). 

This paper is organized as follows. We first introduce our model of the primordial non- 
Gaussianities and discuss the basic statistical properties in Sec. 2. We then compute the 
halo power spectrum following the peak-background split argument in Sec. 3. We discuss 
the two relevant parameters for the feature in the scale-dependent bias and show how they 
are useful to distinguish different models in Sec. 4. We give a brief summary of the paper 
in Sec. 5. Throughout the paper, we adopt the best-fit flat ACDM cosmology to the seven- 
year observations of WMAP [3] and adopt A s = 2 A3 x 10 -9 , n s = 0.963 for the amplitude 
and tilt of the power spectrum of the curvature perturbations normalized at wavenumber 
fco = 0.002 Mpc -1 , O m = 0.265, f2b — 0.0448 for the density parameters of total matter and 
baryon and h = 0.71 for the Hubble parameter. We compute the transfer function from the 
curvature perturbation to the matter density by CAMB [54] using these parameters. 

2 The model 

In this section, we describe our model of the non-Gaussian curvature perturbations and show 
their basic statistical properties. This section gives useful quantities to derive the clustering 
of halos in later sections, such as the cumulants of the density field. We first introduce our 
parameterization for the curvature perturbations in Sec. 2.1, and compute the polyspectra 
(Sec. 2.2) and cumulants (Sec. 2.3). 

2.1 Local-type non-Gaussianities with two fields 

In this subsection, we introduce our models of the non-Gaussian primordial perturbations. 
Throughout the paper, we restrict our attention to the adiabatic density fluctuations, whose 
statistical properties are solely transferred from the primordial curvature perturbations: 



where M(k, z) denotes the matter transfer function, and is normalized such that 5k(z) gives 
the linear overdensity field at redshift z. We hereafter omit the redshift dependence, but 
it always comes through M(k,z). Then, we consider the local- type non-Gaussianities in 
the curvature perturbations originating from two Gaussian fields. We expand the curvature 
perturbations into a Taylor series: 



where xi anci X2 are two statistically independent auxiliary Gaussian fields, (X1X2) — 0. In 
the above, the indices i and j in the summation run over positive integers with i + j > 2. We 



5k(z) =M(k,z)( k 





(2.2) 
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assume that the linear terms in Eq. (2.2) are dominant contributions and the higher-order 
terms give little corrections to the total curvature perturbations. 

After the main discussion starting with the two-field model in Eq. (2.2) we will show 
how the result can be generalized to multi-field models with N > 2 independent Gaussian 
fields. Our final results for the halo power spectrum do not depend on the number of fields 
as will be shown in Sec. 3.2. 

2.2 Polyspectra 

Let us first compute the polyspectra based on Eq. (2.2). We focus on models in which the 
two fields, Xi an d X2, have the same spectral index, and parameterize their power spectra as 

P xt {k) = (l-a)P f (fc), (2.3) 
P X2 (k) = aP c (fc), (2.4) 

n 2 / u \ "8-1 

Note that Eq. (2.2) reduces to the single-field model (1.1) when a = 0. We can generalize 
the model (2.2) into scale-dependent non-Gaussianities by setting different tilts for the two 
fields, and these models will be discussed elsewhere. 

We then consider the higher-order polyspectra of £ at the tree-level. The model (2.2) 
reproduces the bispectrum of the curvature perturbations expressed by Eq. (1.2), with the 
amplitude parameter f n \ being the following form: 

5 

/ni = g [(1 - «) 2 c 2 ,o + a(l - a)ci,i + a 2 c , 2 ] • (2.6) 

Similarly, the model gives the trispectrum written as the sum of Eqs. (1.3) and (1.4), and 
the relevant parameters read 

25 

= 7T [(1 — a) 3 C3,o + a(l - a) 2 c 2 ,i + a 2 (l - a)c 1>2 + a 3 c 0y3 ] , (2.7) 
Tni = 4(1 - afc\ Q + q(1 - a)cj A + 4a 3 Co 2 + 4a(l - a)c lt i [(1 - a)c 2>0 + ac , 2 ] . (2.8) 

The trispectrum has two terms, one comes from the cubic coupling, Cj with i + j = 3, 
and is characterized by g n \, while the other originates from Cjj with i + j = 2 (i.e., the r n i 
term). If one specifies the values of / n i, g Q \ and r n i, the primordial trispectrum as well as the 
bispectrum are uniquely determined. 

In this paper, we further discuss the primordial non-Gaussianities beyond the trispec- 
trum level to confirm the generality of the conclusion. In order to do so, we compute the 
tetraspectrum (the fifth order spectrum) of the curvature perturbations. In addition to the 
contribution from h n \ in Eq. (1.1), we find two more terms arising from lower-order non- 
Gaussianities (i.e., / n i and g n \). They are shown in Appendix A. As a specific example, we 
investigate the effect of quartic-order non-Gaussianity characterized by h n \, and the results 
of iV-body simulations will also be presented in paper II. 

2.3 Cumulants 

The cumulants of the density field are basic quantities useful in analytically estimating the 
halo number density (see Sec. 3.2 below). In this subsection, we compute the low-order 
cumulants and give their fitting formulae. 
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We begin by denning the smoothed density field in Fourier space: 

Srjc = W R {k)8^ (2.9) 

where W R is the Fourier transform of the window function at scale R, and we adopt the top 
hat function. We then define its n-th order cumulant (or, the connected moment) as 

4 n) = <«a>c (2.io) 

We introduce the windowed transfer function, Ai R (k) = W R (k)M(k), which directly con- 
nects the smoothed density contrast to the curvature perturbation: 

5 R , k = M R (k)Cu- (2-11) 
Using Gaussianity of xi an d X2> the leading order non- vanishing cumulants are computed as 



H, — ,tn>3X> — rfr. 



n 

m=l 



\ 

II [Mr * (xTxt - (xTxt))} ) , (2-12) 



\m=l 



where the product * represents a convolution with the windowed transfer function: 

(•Mr * A)(x) = J d 3 yM R (x - y)A(y), , (2.13) 

with A^_r(x) being the Fourier counterpart of A4 R (k). In Eq. (2.12), the summation runs for 
In non- negative integers, i m and j m , which satisfy 

«m + Jm>l, for m = 1,2, ... ,n, (2.14) 

n 

Y J ( i rn+3rn)=2{n-l). (2.15) 
m=l 

They can be explicitly rewritten by integrals of the polyspectra of £, already computed in 
the previous subsection, after multiplying the window function. For example, the first three 
cumulants, the variance R R , the skewness R^ and the kurtosis R^ , are given by 

4 2) = / -^Mr(p) 2 P ( (p) = 4, (2-16) 

4 3) = / ^^^(p)^i?(Q)^(|p + q|)^c(p>9>lp + ql) = /ni4 3 / nl ' (2-i7) 

_g) = J d3 P^q d3r M R (p)M R (q)M R {r)M R (\p + q + r|)T c (p, q, r, -p - q - r), 

^Snl^L+^L (2.18) 

where f^ Rg l and Rr t l originate from the trispectrum in Eq. (1.3) and (1.4), respectively. We 
find that the following fitting formulae for the skewness and kurtosis accurately reproduce 
the numerical integrals of the above equations for our fiducial cosmology at z = 0: 

kf f -3.20 x 10~ 4 a| 141 , (2.19) 



*R,f, 

4 4 ) 9ni ^5.81xl0- 8 4 297 , (2-20) 
~ 2.06 x 10~ 7 4 288 - (2-21) 
See also [55] and [56] for similar fitting formulae. We also derive analogous expression for 
the fifth-order cumulant, k r . It can be found in Appendix A. 
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3 Analytical prediction of the halo clustering 

So far, we have focused on the basic properties of the underlying density field as well as the 
curvature perturbations. Now, let us consider a simple model of the halos in the presence of 
general local-type non-Gaussianities described by Eq. (2.2). 

We follow the peak-background split formalism, and generalize it to our model, Eq. (2.2) 
in Sec. 3.1. There, we define the conditional cumulants and show how these cumulants are 
affected by primordial non-Gaussianities. After that, we discuss the clustering properties of 
halos based on the calculation of the conditional number density of them in Sec 3.2. 

3.1 The peak- background split in the presence of non-Gaussianities 
3.1.1 long and short mode decomposition 

The peak-background split formalism had been originally introduced to model the biased 
clustering of peaks, halos, or galaxies which are associated with the initial peak regions in 
the Lagrangian space [57-60]. This formalism gives a simple framework to compute the bias 
of the collapsed objects with mass M by considering the conditional probability that a density 
contrast with a mode of scale R exceeds a critical density 5 C under the background density 
field modulated by long-wavelength modes. Here, the mass and the scale are related as 

where p denotes the cosmic mean density. Effectively, the presence of long modes can be 
interpreted as a local shift in the threshold density: 

5 C S c (q) = 5 C - fc(q), (3.2) 

where the coordinate q indicates that we are working in the Lagrangian space. We adopt 
6 C = 1.686 following the result of the spherical collapse model, and leave the discussion on 
its more accurate modeling to paper II. Using Eq. (3.2), one can relate the long mode to the 
halo number density at that position with a help of the halo mass function. 

In the presence of primordial non-Gaussianities, Eq. (3.2) is not the only effect from 
the long modes, since Fourier modes are not independent of each other from the beginning. 
We will show how this mode coupling results in the statistical properties of the short modes 
that are responsible for the formation of halos. 

We start with decomposition of the two auxiliary Gaussian fields into sums of the long 
and short wavelength contributions: 

Xi(q) = XM(q) + Xi,«(q)> * = 1, 2. (3.3) 

Then, keeping only the linear modulation from Xiti the short mode for the density fluctuation 
can be computed as 

<5 s (q) = £ [aj + Acij(q)] (Mr * \,„WKq). (3.4) 

y 

Acij(q) = (i + l)c i+ ij Xl,<(<l) + (j + l)cy+i X2/(q), (3.5) 

where ci o = co ; i = 1 and the convolution * is defined by Eq. (2.13). Interestingly, one can 
see that a (n + l)th-order coupling constant in £ (i.e., Cj with i + j = n + 1) modifies 5 S at 
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the nth order. For instance, in case of quadratic non-Gaussianities characterized by f n \, the 
short mode is affected by the long mode at the linear order (i.e., the amplitude of the power 
spectrum is modulated), and this gives the scale-dependent bias [35]. We will see how the 
statistical properties of the conditional density field (3.4) result in the characteristic feature 
in the halo clustering shortly. 

3.1.2 conditional cumulants 

In Sec. 2.3, we have computed the cumulants of the density field in the presence of primordial 

-(n) 

non-Gaussianities. The cumulants, denoted by k r , are given as the volume averages of the 
products of the density contrast over the entire universe. One of the key idea in the peak- 
background split formalism is to consider the conditional cumulants of the short mode at an 
arbitrary position q in the presence of the long modes. These cumulants are defined as the 
ensemble average of 5^ while keeping the long mode at that Lagrangian position (e.g., Xi/il) 
in our case) fixed: 

4 n) (q) = (<^W (3-6) 

Then, we can easily show how the long mode, which couples with the short mode due to the 
primordial non-Gaussianities, modify the cumulants of the short mode locally. 

In Eq. (3.4), we can see that the nonlinear coupling constant in the short mode of the 
density contrast can be obtained by replacements of the coupling constants to the effective 
ones 

c ij -> c *j(q) = c t,j + Ac i,j(q)> ( 3 - 7 ) 

which vary from position to position proportional to a certain combination of Xli an d X2£- 
We can obtain the conditional cumulants by applying the above replacements to Eq. (2.12) 
derived for the unconditional cumulants. Let us consider first the conditional variance of 
the short mode fluctuations of the density contrast. We define the local rms fluctuation, 
0"ij(q) = [fi^(q)] 1 / 2 , and this is given by 

o-fl(q) = {i + [ 2 (! - ") c 2,o + aci,i]xi,*(q) + [(1 - «)ci,i + 2qc , 2 ] X2,e(<l)jv R , (3.8) 

where we keep only linear corrections from xi,e an d X2£- Similarly, we can compute the 
higher-order conditional cumulants of the short mode when the long modes are fixed. Again, 
we keep the leading-order corrections arising from Xi,£i an d obtain 

4 3) (q) = [u + A/„i(q)]/4 3 L' ( 3 - 9 ) 

4 4) (q) = [9nl + A 5n i(q)]«g ) flni + [r nl + Ar nl (q)]4 4) Tni , (3.10) 

for the skewness and kurtosis. Here we define the effective shifts in the coupling constants, 
A/ n i(q), A<7 n i(q) and Ar n i(q). They can be explicitly written down as linear combinations 
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of Xi^(q) and X2,e( c l) : 



A/m(q) = g{ [C 1 - «) 2 ( 4 4o + 3c 3,o) + «(1 - a)(2ci,ic 2l o + <? x ,\ + ^2,1) 



Xi/(q) + 



1 o 2 



Agni(q) 



5 
3~ 

+a 2 (2c , 2 ci,i + ci )2 ) 
25 r r 

— | 2(1 - a) 3 (3c 2 ,oC3,o + 2c 4j0 ) + a(l - a) 2 (4c 2 ,oc 2 ,i + 01,102,1 + 3c 3 ,i) 
+2a 2 (l - a) (02,0^,2 + 01,101,2 + 02,2) + a 3 (3ci,ic ,3 + 01,3) 



X2,e( c 0fi 



(3.11) 



+ 



1 o 2 



Ar n i(q) 



X2,i 



(q)}, 



(3.12) 



(l-a) d (24 + 3c2,oc 3 ,o) 
+4q(1 - a) 2 (2c 2 , ci,i + 2c2, ci,i + ci,ic 2 ,i + 3ci,ic 3 , + 2c 2 ,oc 2 ,i) 
+2a 2 (l - a)(c? 1 + 4c 2 ,oCi,ic ,2 + 2c 2 iC , 2 + 2ci,ic 2 ,i + 4co, 2 c 2 ,i + 2ci,iCi, 2 



+8a 3 (ci,iCQ2 + co, 2 ci, 2 ) 



Xii(q) + 



1 o 2 



X2J? 



;q)}, 



(3.13) 



where the coefficients for X2/( x )> denoted as [1 -H- 2], can be obtained by performing permu- 
tations c m , n c n , m and a -H- (1 — a) to the corresponding coefficients of xi/( x )- See also 
Appendix B for the results for single-field non-Gaussianities (1.1). 

An important point here is that the nth-order correction in Eq. (2.2) modifies the nth- 
order conditional cumulant, while the unconditional cumulant at the same order is not af- 
fected. For instance, if one is interested in quadratic models with f n \ 7^ 0, the first non-trivial 
contribution to the conditional cumulants appears in the variance, while the unconditional 
variance does not depend on f n \. Similarly, when g n \ 7^ 0, the conditional skewness has a 
linear dependence on the parameter g n i (see also [30, 43]), while this parameter gives the 
amplitude of the unconditional kurtosis. 

3.2 Halo number density and clustering 

The halo mass function is known to be almost universal and solely determined by the pa- 
rameter v = 5 c /aR, where 5 C denotes the critical density above which halos form at z = 
[61]. The mass function is usually characterized by f(M,z): 



n(M, z) 



dM 



f(M,z) 



P 
M 2 



d In Up} 



dlnM 



(3.14) 



where we denote by n(M, z) the number density of halos with mass M at redshift z, and 
p is the cosmic mean density [62]. When the halo mass function is universal, the function 
/(M, z) depends only on the peak height v. However, as we are interested in corrections 
from the higher-order cumulants, we slightly relax the assumption of the universality. We 
allow a dependence on the density cumulants at any order and write 



f(M,z) =fla R ,4\R < £ ) 



(3.15) 



Here, the argument R has a one-to-one correspondence to the halo mass M via Eq. (3.1). 

We then assume that the halo number density field is a local process in Lagrangian space: 
given the conditional cumulants defined in the previous subsection, we can deterministically 
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predict the number of halos at that (Lagrangian) position. More specifically, we assume 
that the halo number density is a function of the conditional cumulants smoothed with 
the physical scale corresponding to the halo mass in the Lagrangian space. We introduce 
/(q|M, z) analogously to its unconditional counterpart, /, in Eq. (3.14), and assume that 
this function can be obtained by replacements of the arguments in / as 

*£> ->«£>(<!), (3.16) 

in addition to the change in the critical density for the formation of halos in Eq. (3.2). We 
have 

f(q\M, z) = f (a R (q), (q), (q), fc(q)) . (3.17) 

This function gives us a simple route to the halo power spectrum. The halo density 
contrast in the Lagrangian space can be calculated as 

XL, \ f{<l\ M ,z) n fa ion 

du (q \M, z) = r 1. (3.18) 

By expanding this equation into Taylor series of Sg, xi£ an d Xit-, an d collecting all the 
leading-order contributions, we obtain the linear halo overdensity field in the Lagrangian 
space. We then map this field to the Eulerian space. At leading order, we have 

8 h (x) = b s 5(x) + b xi xi (x) + b X2 X2(x), (3.19) 

with coefficients defined through the derivatives of the mass function: 

(3.20) 



91nf din an ^ din f 

t »-55^ftsr + Ss^*S7' s ^' 2 ' (3 ' 2I) 



In Eq. (3.19), we dropped the subscript i for notational convenience. Note that the first term 
in Eq. (3.20) comes from the mapping to the Eulerian space. Note also that the derivatives 
of the conditional cumulants with respect to Xit can De written down explicitly in terms of 
Cij and a by using the expressions derived in Sec. 3.1.2. Taking the average of the square of 
Eq. (3.19) in Fourier space, we finally arrive at the halo power spectrum: 

P h (k) = b 2 & P s (k) + 2r UF h\P sc (k) + blP c (k), (3.22) 

where the density and density-curvature power spectra are defined as 

P 5 (k) = M 2 (k)P ( (k), P sc (k) = M(k)P c (k), (3.23) 

and we parameterize the bias coefficients by 



& c = ^(l-a)&2 i+a &2 2 , (3.24) 
_ (1 - a)b xi + ab X2 

TMF = , =■ (3.25) 
/(l_ a)6 2 i+a6 2 2 
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As we will discuss later in more detail, the parameter characterizes the amplitude of 
the non-Gaussian corrections to the halo power spectrum. On the other hand, the other 
parameter tmf captures the shape of the scale dependent bias, and plays a key role to 
distinguish the models for the primordial non-Gaussianities based on multi fields from those 
based on a single field. 

It is interesting to see that a wide class of models described by Eq. (2.2) end up with 
the same result for the halo power spectrum with only two parameters responsible for non- 
Gaussian corrections (b^ and tmf) in addition to the usual (Gaussian) bias factor, 65. More- 
over, our result can be generalized to models with more than two independent fields very 
easily. If we employ N statistically independent Gaussian fields, Xm an d expand the curva- 
ture perturbations as 

C(x) = J^Xn(x)+ 

the halo density contrast of Eq. (3.19) becomes 

*h(x) = M(x) + ^6 XnX n(x), (3.27) 

n 

with the index n running from unity to N. From this equation, it is clear that Eq. (3.22) 
still holds as long as all the N fields have the same spectral index. In this case, the bias 
coefficients are given by 

k = ./£«n&2„> ( 3 - 28 ) 



nviF = 




(3.29) 



with the fraction in the power spectrum amplitude of the n-th field denoted as a n that satisfy 

X)«n = 1 - ( 3 ' 3 °) 

n 

Equation (3.27) can be regarded as a generalization of the multi-variate biasing in case of 
the single- field models discussed in [28]. 

We show an example halo power spectrum in Fig. 1. We can see how each of the three 
terms in Eq. (3.22) contributes to the total power spectrum. The non-Gaussian corrections 
are prominent on large scales ( £j 0.01/iMpc" 1 ). The parameters, and tmf as well as b$ 
can be determined from observed power spectrum with a sufficient survey volume thanks to 
the different ^-dependence of the three terms. We propose to use these parameters as direct 
observables from the measurements of the scale-dependent halo bias instead of popular non- 
Gaussian parameters such as f n \ and g n \, because the latter set of parameters have strong 
degeneracy in many cases. We will see how our new parameterization helps us to test different 
class of models for the primordial non-Gaussianities in the next section. 



IlXn(x) 




(3.26) 
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k [h Mpc" 1 ] 

Figure 1. Example plot of the halo power spectrum (Eq. 3.22) in case of b$ = 2, = 100 and 
fMF = 1. We plot the total power spectrum by solid line, while the three terms in that equation are 
respectively shown in dashed, dotted and dot-dashed lines. 




3.3 Recovery of previous studies 

It is worth comparing our results with the predictions in the literature. In the presence of 
the primordial non-Gaussianities based on a single field, it is straightforward to show that 
rMF = 1 or — 1 (see Sec. 4.1 for more discussion on tmf)- Then, Eq. (3.22) can be rewritten 

as 

P h (k) = [b 5 + Ab(k)] 2 P s (k), (3.31) 

Ab(k) = sign^MF)^-^" 1 ^)' ( 3 - 32 ) 

where we denote the sign of A as sign(^4). In this case, one can obtain the scale-dependent 
bias in the halo power spectrum simply by replacing a constant bias factor with a fc-dependent 
one: b$ — > b$ + Ab(k). All the effects from the primordial non-Gaussianities are encoded in 
b^, the amplitude of the scale-dependent bias. This parameter in case of the single-field 
primordial non-Gaussianities (1.1) is explicitly given in Appendix B. 

Let us further restrict the model to the quadratic non-Gaussianities (i.e., Ad-type). We 
find that the first term in Eq. (3.21) is the dominant contribution to b^, and the derivatives 

in) 

of cumulants k r with n > 3 are small. By dropping these small corrections, and using that 
the derivative with respect to 5 C in Eq. (3.20) can be converted to that with respect to <tr 
under the assumption of the universal halo mass function, we recover 

Ab(k) = ^US c (b s - l)M-\k), (3.33) 
5 

derived by [23]. The correction terms from higher-order cumulants might be important if 
one needs a very accurate modeling for the scale-dependent bias. However, since these terms 
strongly depends on the non-Gaussian mass function and we adopted a rather simplified one 
described in Appendix C in this paper, we left further discussions on the accurate modeling 
in paper II. 



We next consider cubic non-Gaussianities described by g n \. In this case, the first term 
of Eq. (3.21) is zero, and thus the dominant effect comes from the derivative of the mass 
function with respect to the skewness. Using the expression for the conditional skewness 
given in Appendix B, and dropping the corrections from higher-order cumulants, we have 



This is identical to Ref. [43]. One can recover the formula for the scale-dependent bias 
derived in [30, 31] by substituting the mass function based on the Edgeworth expansion 
to this expression (see Appendix C). Indeed, Ref. [43] demonstrated that the A-body data 
can be accurately explained by evaluating the derivative of the mass function directly from 
simulations. 

Let us put some comments on non-Gaussianities based on multi-field initial conditions. 
Refs. [46, 47] compute the clustering of halos in a two-field inflationary model using the peak- 
background split. They employ two fields, one is Gaussian and the other is non-Gaussian 
with a quadratic correction (i.e., / n i-type), and thus our model (2.2) includes their model. 
However, their main focus is on the stochasticity between the halo and the underlying matter 
density fields. Although our derivation in this paper is equivalent to theirs, our focus is more 
on the auto-power spectrum of halos. It is straightforward to show that the stochasticity 
arises between the halo and the matter fields in our models when a / regardless of the 
order of the non-Gaussian corrections, and this is another unique feature to test multi-field 
non-Gaussianities. 

Note also that our analytical calculation of the halo power spectrum is valid only at the 
linear order. Although we are interested in the clustering on large scales (say, ^ l/i~ 1 Gpc) 
and linear theory is expected to work well, any nonlinearity (e.g., nonlinear gravitational 
growth or nonlinear bias) can in principle modify our results. For example, the nonlinearity 
in the halo bias can source a similar scale dependence as shown in [24] (and see also [30] for 
a comparison among different derivations of the scale-dependent bias including the effect of 
nonlinear bias in the thresholding scheme). As another example, the matter power spectrum 
is affected by a loop correction in case of g n i-type models, which does not vanish at large 
scale limit [42]. Another important assumption in our analysis is the locality of the formation 
of halos in Lagrangian space. These issues will be discussed in paper II with fully nonlinear 
treatment based on cosmological A-body simulations. 

4 Testing non-Gaussian models with observation 

We have shown that our model (3.26) generally results in the halo power spectrum Eq. (3.22). 
The final result of the halo power spectrum contains only three parameters, b$, and rMF) 
and they can be determined or constrained directly from the observation using the different k- 
dependence of the three terms in Eq. (3.22). Following this result, we discuss what properties 
of the primordial non-Gaussianities can be tested by measurements of these parameters. 

4.1 Cauchy-Schwartz inequality and multi-field non-Gaussianities 

We discuss what we can learn from the parameter tmf defined in Eq. (3.29). By simple 
calculation, one can show that tmf satisfies the following inequality: 



Ab(k) = -g nl k 



K fm,R 



d In/ 



M~ x {k) 



(3.34) 



ryre\ < 1 
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where the equality holds when the curvature perturbations originate from a single field. In 
multi-field models, the above equality holds only when all the coefficients b Xi are the same. 
Note that this inequality is the Cauchy-Schwartz inequality among Gaussian fields, Xi- 

Measurement of this parameter opens a new window to distinguish between single-field 
and multi-field primordial non-Gaussianities. If the absolute value of tmf measured from 
the observations of the galaxy clustering significantly deviates from unity, we can rule out 
single- field models of the primordial fluctuations. Interestingly, this inequality is shown to 
be identical to the Suyama-Yamaguchi inequality in case of the quadratic non-Gaussianities 
under some approximations. Let us focus on two-field models (2.2) for simplicity. As we 
have already mentioned, the dominant contribution to the bias coefficients in Eq. (3.21) is 
the first term in the presence of quadratic non-Gaussianities, Cjj with i + j = 2. If we drop 
the higher-order corrections, we have 

b Xl « [2(1 - q)c 2i0 + aci,i] ttt^-, (4.2) 

dlna R 

d In/ 

Substituting these expressions into Eq. (3.25), we finally have 



r M F » i~F=> ( 4 - 4 ) 



with a help of the definitions of / n i and T n \ in Eqs. (2.6) and (2.8). Note that Eq. (4.4) still 
holds in more general case of Eq. (3.26). Thus, we can directly test the Suyama-Yamaguchi 
inequality from the observed power spectrum. Note that our parameter tmf is equivalent to 

— 1/2 

A NL introduced in [12] under this approximation. Since our inequality is valid for models 
with higher-order non-Gaussianities, it can be regarded as a generalization of the Suyama- 
Yamaguchi inequality, and can be directly tested from the observed halo power spectrum. 

The parameter tmf controls the shape of the scale dependent bias. This is illustrated in 
Fig. 2. In the left panel, we show how the halo power spectrum depends on tmf when bg and 
are fixed. The solid lines correspond to models of primordial non-Gaussianities originating 
from a single-field with f n \ = ±100, while the dashed lines show the halo power spectra from 
multi-field non-Gaussian initial conditions. We also plot by symbols the expected error bars 
from a 100^ 3 Gpc 3 survey at z = 1. We employ the formula by [63] for the \-a error bars 
that is valid for Gaussian fields: 



[AP g (kf 



-^modc 



(4.5) 



where iV moc j e denotes the number of independent Fourier modes. The fiducial model shown 
as the error bars in the left panel corresponds to f n \ = 20 and r n i = 14, 400, if we assume 
that the non-Gaussianities are quadratic. The value of f n \ is small enough to satisfy the 
current observational bounds, while that of r n \ is much larger than (36/25)/^, which might 
be possible to detect from future observations. 

From the figure, we can see the largest difference among the five lines at k ~ 0.005/iMpc -1 , 
where the Pg^ term in Eq. (3.22) has the largest impact. On the other hand, the they con- 
verge at the small-fc limit since we adopt the same value of for them. We need to measure 
the broadband shape of the power spectrum covering this entire scale at £ h~ 1 Gpc to obtain 
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Figure 2. Dependence on the shape parameter tmf of the halo power spectrum. Left: analytical 
results of the tmf dependence at z = 1. Lines show models with different tmf while b$ = 2.1 and 
£>C = 210 are fixed. Expected error bars on the galaxy power spectrum from a 100/i~ 3 Gpc 3 class 
survey are also shown assuming tmf = 0.2 as the fiducial value. We assume that the typical host halo 
mass is 10 13 h~ 1 M Q (this gives the value of bs shown in the figure), and the number density of halos 
is Hh = 10~ 4 /r 3 Mpc™ 3 . Right: examples for halo power spectra measured from iV-body simulations. 
We fit the data by model (3.22) to find the best-fit parameters, bs, b^ and tmf- They are shown in 
the figure for each model. 



meaningful constraints on tmf- The planned surveys such as EUCLID [64] can measure the 
galaxy power spectrum quite accurately even at such a large scale thanks to their huge survey 
volumes. They will provide great opportunities for this kind of tests. 

We also show in the right panel of Fig. 2 the measured halo power spectra from N- 
body simulations at z = 1. We plot the results obtained in simulations starting with two 
different non-Gaussian initial conditions; a single-field model (circles), and a two-field model 
(triangles). The input parameters of these simulations are C2,o = 60 and a = for the single- 
field model, and C2,o = 120, co,2 = —120 and a = 1/2 for the two-field model. The derived 
non-Gaussian parameters for the two simulations are (/ n i, r n i) = (100, 14400) and (0, 14400), 
respectively. Thus, the two simulations are expected to have roughly the same values for bg 
and &£. We fit the measured halo power spectrum by Eq. (3.22) to obtain bs, and fmf, 
which are shown in the panel. The best-fit value of tmf hi case of the single-field model 
shown in the figure is consistent with unity within the statistical uncertainty, as expected. 
On the other hand, the multi- field model clearly disfavors |tmf| = 1- These results basically 
confirm our analytical expectations. See paper II for more detail on the simulations and 
analyses. 

It is interesting to note that in Ref . [45] , the authors reach a similar conclusion for testing 
/ n i and r n i through the shape of the halo power spectrum but from a different derivation based 
on the clustering of peaks above a threshold. Our results agree with theirs qualitatively. One 
advantage of our formulation is that our inequality (4.1) is always valid with non-Gaussian 
corrections at an arbitrary order and can be used to test multi-field models for the primordial 
fluctuations in general. 



- 15 - 




1 2 3 4 51 2 3 4 5 



b§ bg 

Figure 3. Dependence on the Gaussian bias factor, bs, of the non-Gaussian correction to the halo 
bias. We plot the combination A n Q = b^/(bs — 1) as a function of bs at redshifts 0.35, 1 and 2. The 
error bars show the 1-a uncertainty obtained by fitting to TV-body simulations, while the lines depict 
our theoretical predictions at z = 0.35 (solid), z = 1 (dashed) and z = 2 (dotted). The left (right) 
panel plots the results for quadratic (cubic) non-Gaussianities characterized by f n \ (g n \). 



4.2 Mass dependence and the order of the non-Gaussian corrections 

We finally discuss how we should distinguish the quadratic non-Gaussianities (i.e., Ad-type) 
from higher-order models. As we have shown in the previous subsection, the shape of the 
halo (or galaxy) power spectrum on large scales only gives the information about the number 
of independent fields through the parameter tmf- However, the / n i-type non-Gaussianities 
are distinguishable from higher-order ones if one examines more than one galaxy populations 
with different bias parameters, in principle. In this section, we discuss how the ratio 

can be used for this test. 

For simplicity, let us focus on the primordial curvature perturbations originating from a 
single Gaussian field (1.1). Figure 3 shows how the amplitude of the non-Gaussian correction 
to the halo bias scales with the Gaussian bias factor for halo catalogs with different masses 
at different redshifts. We plot A u q measured from iV-body simulations (error bars) and 
estimated by our analytic calculation (lines). We adopt two models for the primordial non- 
Gaussianities, and show the results for f n \ = 100 in the left panel and g n \ = 5 x 10 5 in the 
right panel. 

As expected from Eq. (3.33) derived in the literature, this ratio remains almost con- 
stant over the Gaussian bias factor (or the halo mass) and redshift in case of quadratic 
non-Gaussianities. In contrast, the scaling of A n Q with bs observed in the right panel is very 
different. A clear increasing tendency is visible from both simulations (symbols) and analyt- 
ical predictions (lines). Interestingly, the redshift dependence of the relation between A u q 
and bs is very weak similarly to the left panel for the / n i-type primordial non-Gaussianities. 
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Since we need an accurate modeling of the derivatives of the mass function with respect to the 
cumulants to compute the analytical curves, it is generally difficult to reach a quantitative 
agreement with simulations. See Appendix C for more detail of the mass function we adopt 
in computing the analytical predictions. 

Although our analytical calculation has a slight offset to the iV-body data, the qualita- 
tive agreement with the A-body simulations is rather encouraging. The difference between 
the two models shown in Fig. 3 originates from the different origins of the scale-dependent 
bias. In the presence of a non-zero / n i, most of the effect comes from the conditional variance 
as was shown in the previous section. In this case, roughly scales as oc (b$ — 1), and 
thus A n Q does not depend on the halo mass. In contrast, we find that the origin of the 
scale-dependent bias is in the conditional skewness in the presence of non-zero g n \. 

Using this feature, we propose the following test. Measure the power spectrum from 
multiple tracers, such as blue and red galaxies, massive and less massive halos at any redshifts. 
Then, fit by Eq. (3.22) to obtain b^ as well as b$ and compare the ratio A n Q measured from 
each tracer. If the two (or more) values of A q q are inconsistent with statistical significance, 
that might be a sign of the primordial non-Gaussianities which originate from higher-order 
coupling (i.e., g n \ or even higher). We should have a constant A n Q independent of the 
properties of the tracers in case of the quadratic non-Gaussianities characterized by f n \. We 
also examine the parameter A u q in the presence of quartic-order non-Gaussianities described 
by h n \ (see Eq. 1.1). We find that the bs dependence of this parameter is similar to what is 
found in the presence of cubic non-Gaussianities. Thus, we conclude that the constancy of 
A u q is generally a key feature for discriminating / n i-type non-Gaussianities from higher-order 
models. 

5 Summary and Discussions 

We computed the halo power spectrum in a class of local-type non-Gaussianities described 
by Eq. (3.26) based on the peak-background split formalism. We found that the resultant 
halo power spectrum generally exhibits a scale dependence in the large scale bias factor 
and is parameterized by three parameters; the ordinary Gaussian bias factor, bs, the shape 
parameter of the scale dependence, tmFj and the amplitude of the scale dependence, b^. 
We showed how one can distinguish between non-Gaussianities based on a single field and 
multi fields through the measurement of tmf • This allows us a direct test of the generalized 
Suyama-Yamaguchi inequality (4.1) from observational data. Another parameter, b^ might 
be useful to distinguish / n i-type models from g n i-type or models with even higher-order non- 
Gaussianities. The parameter ^4 n G, a combination of b^ and b§, is shown to be almost constant 
in case of / n i-type non-Gaussianities, though it can have a strong dependence on the mass 
of halos in general. Thus, we propose to measure this parameter from multiple tracers to 
search for a possible dependence of A n Q on the halo mass to test higher-order non-Gaussian 
models. 

There are several concerns on the systematic error and possible degeneracy with other 
models for the primordial non-Gaussianities. Since our analytic model is based on the leading- 
order calculation, any loop-correction can be a source of systematic error. This can be a 
problem for both of our two tests, if the correction is large. We thus have to test our model 
against simulations. A detailed comparison with iV-body simulations will be presented in a 
separate paper. 
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Another systematic effect, though it is avoidable by a careful theoretical modeling, is 
induced from a mass dependence in the critical density for the halo formation, so called the 
moving barrier [65]. According to the moving barrier argument, the critical density 5 C might 
not be a constant, and should be replaced by 5 C (M) with an explicit halo-mass dependence. 
In the presence of non-zero f n \, the ratio A n Q can be approximated as 



when we neglect terms from higher-order cumulants. This equation implies that the mass 
dependence in 8 C (M) can be misinterpreted as a sign of higher-order non-Gaussianities. We 
will discuss in paper II how this matters using iV-body simulations, and show that in the 
presence of non-zero f n \, we have no clear evidence of the scale dependence in 5 C (M). 

The halo assembly bias is another source of systematics. Ref. [66] showed that the 
amplitude of the scale-dependent bias for a given halo mass can be different when they divide 
the halo sample into subsamples with different halo assembly histories. This can be a problem 
when the observed galaxies form in halos with a certain merging history preferentially. The 
problem of non-Gaussian halo assembly bias may also be interpreted as a dependence on the 
assembly history of the barrier function S C (M). We need a careful modeling of the barrier 
function as a function of the assembly history in addition to the mass. 

Finally, scale dependent f n \ models can also result in a mass dependent ratio A q q. In 
Ref. [49], the authors discuss the scale-dependent bias in models with k dependence in the 
coupling parameter f n \(k). They show that the magnitude of the non-Gaussian correction to 
the bias in the power spectra of massive and less massive halos depends on f n \(k) evaluated at 
wavenumber roughly corresponding to the mass scale of the halos (see also [30, 31, 48]). Thus, 
this dependence can mimic the mass dependence of ^4 n c that we found in the presence of 
higher-order non-Gaussianities. However, by detecting a mass dependence of this parameter, 
we might be able to falsify the simplest non-Gaussian model with a constant / n i, and thus 
this test is still interesting. 

This work is the first step towards establishing a methodology to test non-standard 
primordial non-Gaussianities beyond the so-called / n i type in the era of 100/i~ 3 Gpc 3 class 
full sky surveys. We expect that our tests for the multi-field non-Gaussianities as well as the 
higher-order ones will be possible with future surveys such as EUCLID. 
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A Fifth-order statistics 

In this appendix, we show the relevant formulae for fifth-order statistics in model (3.26). In 
actual calculation of the bias coefficients (3.21) in the text, we truncate the series at the 
fifth-order cumulants. The numerical results shown in Fig. 3 are based on the formulae in 
this section. 

First, the tetraspectrum of the curvature perturbations reads 

T f (5) (ki,... ,k 5 ) = ^/ inl [P c (fc 1 )P c (A; 2 )P c (A ; 3)P c (A ; 4) + (4perms.)] 

+Tni ] [Pd k i) P d k 2) P d k 3)Pd\ k i + k 4l) + (59 perms.)] 

+r^ U) [P ? (A; 1 )P c (A; 2 )P c (|k 1 + k 3 |)P ? (|k 2 + k 4 |) + (59perms.)] .(A.l) 

In the above, the coefficients, and t^j are respectively given as 

(II) _ 324 (in) _ 216 3 , 

when Q originates from a single field as in Eq. (1.1), and they can be different in multi-field 
cases. Corresponding to the three terms for the tetraspectrum, we give the fitting formulae 
for the fifth-order cumulant of the smoothed density field at z = 0: 

(A.3) 

(A.4) 
(A.5) 

(A.6) 

The fifth-order conditional cumulants are summarized in the next appendix. 
B Exact expressions for single-field non-Gaussianities 

We mainly focus on non-Gaussian curvature perturbations based on multi fields in the text. 
It might be convenient to show the relevant formulae for the single field model (1.1). 

We start with the conditional cumulants. First, substituting a = and with a help of 
Eq. (2.6), the conditional variance of the short-mode fluctuations in Eq. (3.8) is obtained as: 
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°R, (B.l) 



where we use the fact that £ = %i at linear order. This is equivalent to the results of earlier 
works in case of / n i-type non-Gaussianities (e.g., [28, 35]). We then write down the shift in 
the effective nonlinear parameters Eqs. (3.11), (3.12) and (3.13). Substituting a = 0, they 
are simplified as 

A/ni(q) = I (4/ n 2 ! + 3<7ni) CG^(q), (B.2) 

Affnl(q) = g (3/nl&l + 2/i n l) ( G ,M) ) ( B - 3 ) 

216. /„, 2 



Ar nl (q) = — /ni(2/ n 2 1 + 35nl)CG,Kq)- (B.4) 
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(5) 

We can obtain similar expressions for the fifth-order cumulant, n R . This is given as 
4 5) (q) = [h nl + AMq)] 4L + + AT^(q 

where 



k (5) + 

Rr (II) + 
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n > T nl 



A/i n i(q) = - ( 8/ n i/i n i + 5i n i ) Cg,£ (q) , 



972 
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^T(q) = ^|/ n 2 , 



A 7 "™ = ^ ( 4 /nl^nl + 3^! + 6/^ ) C G ,*(q), 



(B.6) 
(B.7) 
(B.8) 



Given the formulae for the conditional cumulants up to the fifth order, it is straightfor- 
ward to derive the bias coefficient, b^. It is 
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where the mass function is give in the next appendix. 



C Non-Gaussian halo mass function 



We need knowledge of the mass function to evaluate the bias coefficients. In this paper, we 
assume a simple mass function and that is described in this appendix. 

We assume that the primordial non-Gaussianities are small, and expand the non- 
Gaussian probability distribution function around normal distribution to obtain Edgeworth 
series as in [67]: 
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where H n (y) is the n-th order Hermite polynomial, and we introduce normalized cumulants, 

= Since we are interested in the primordial non-Gaussianities arising from the 

coupling in the curvature perturbations higher than the quadratic order characterized by / n i, 

_(5) 

we take account of the cumulants up to the fifth, k r , such that we can discuss the effect 
from the quartic correction, h n \. We thus keep the series up to the fourth order in terms 
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of the Edgeworth expansion (see Eq. C.l). We follow the argument by [67] and obtain the 
non-Gaussian halo mass function relevant to this order: 



f(M,z) PBF nG (u) I dcf l 1 dcg } 1 
— r — — — z ti2\y) H 
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# 4 (f) 
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144 dlnafl i/ 61 J 1296 dlna fl i/ 



1 d(c^yi 



(C.2) 



where /g denotes the Gaussian mass function. Note that as long as we consider the local-type 

— ( n ) 

non-Gaussianities, the derivative terms of C R are usually very small given their weak ur 
dependence (see Eqs. 2.19, 2.20, 2.21, and also Eqs. A. 4, A. 5, A. 6), and the first term gives 
the dominant correction. Note also that these derivative terms result in the correction terms 
in the bias parameters derived by [30, 31]. 

Although the above ratio is derived based on the Press-Schechter formalism [62], we 
simply replace the Gaussian mass function with one recently calibrated by a large set of 
iV-body simulations by [68]: 



/ G (M, z) = fmCE(o-R, z) = A{z) 



a 



-a(z) 



+ b(z) 





r <*)] 


exp 









(C3) 



where A(z) = 0.58(1 + z)~ - 13 , a(z) = 1.37(1 + z) -0 ' 15 , b(z) = 0.3(1 + z)- - 084 an d c(z) = 
1.036(1 + z) -0 ' 024 . This mass function is shown to agree with iV-body simulations over five 
orders of magnitude in the halo mass. 

Although the mass function we adopt in this paper to incorporate the effect of the 
primordial non-Gaussianities might be too simplistic, one can easily refine the final result 
by replacing the mass function to a more accurate one. See e.g., [29, 32, 33, 42, 69-71] for 
various elaborated halo mass functions for non-Gaussian initial conditions. 
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